* Create results for Table A5 (col 1 first then col 2)
clear all


* First split doctors in final sample into quartiles based on volumes across the whole period
use "$savedata/masterdata.dta", replace
keep if sample25==1
gen vol = vol25

collapse (sum) one, by(doctor_id)

xtile vol_n4 = one, n(4)

save "$savedata/doctor_vol.dta", replace

* Then run analysis
use "$savedata/masterdata.dta", replace
keep if sample25==1
gen vol = vol25

cap drop _merge
merge m:1 doctor_id using "$savedata/doctor_vol.dta"

preserve

reghdfe survive30 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

foreach x in 1 2 3 4{
xtreg residual c.std_ami3 c.std_nonami3 if vol_n4==`x', fe i(doctor_id)
predict docfe`x', u
gen sigma_u`x' = e(sigma_u)
gen sigma_e`x' = e(sigma_e)
}

collapse (mean) docfe* residual sigma* vol, by(doctor_id)

xtile vol_n4 = vol, n(4)
egen mean_docfe = mean(docfe), by(vol_n4)

gen sigma_u_n4 = sigma_u1 if vol_n4==1
gen sigma_e_n4 = sigma_e1 if vol_n4==1
foreach x in 2 3 4{
	replace sigma_u_n4 = sigma_u`x' if vol_n4==`x'
	replace sigma_e_n4 = sigma_e`x' if vol_n4==`x'
}

gen signal_2step_n4 = ((sigma_u_n4^2) / ((sigma_u_n4^2) + ((sigma_e_n4^2)/vol)))

gen adj_docfe = docfe*signal_2step_n4 + (1-signal_2step_n4)*mean_docfe

su docfe
matrix observe_std = r(sd)
su adj_docfe
matrix observe_adjstd = r(sd)

su docfe
matrix observe_var = r(Var)
su adj_docfe
matrix observe_adjvar = r(Var)

su docfe, det
matrix observe_p10 = r(p10)
su adj_docfe, det
matrix observe_adjp10 = r(p10)

su docfe, det
matrix observe_p25 = r(p25)
su adj_docfe, det
matrix observe_adjp25 = r(p25)

su docfe, det
matrix observe_p50 = r(p50)
su adj_docfe, det
matrix observe_adjp50 = r(p50)

su docfe, det
matrix observe_p75 = r(p75)
su adj_docfe, det
matrix observe_adjp75 = r(p75)

su docfe, det
matrix observe_p90 = r(p90)
su adj_docfe, det
matrix observe_adjp90 = r(p90)

restore

capture program drop myboot2

program define myboot2, rclass

* draw my bootstrap sample
preserve
bsample, strata(doctor_id)

reghdfe survive30 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

foreach x in 1 2 3 4{
xtreg residual c.std_ami3 c.std_nonami3 if vol_n4==`x', fe i(doctor_id)
predict docfe`x', u
gen sigma_u`x' = e(sigma_u)
gen sigma_e`x' = e(sigma_e)
}

collapse (mean) docfe* residual sigma* vol, by(doctor_id)

xtile vol_n4 = vol, n(4)
egen mean_docfe = mean(docfe), by(vol_n4)

gen sigma_u_n4 = sigma_u1 if vol_n4==1
gen sigma_e_n4 = sigma_e1 if vol_n4==1
foreach x in 2 3 4{
	replace sigma_u_n4 = sigma_u`x' if vol_n4==`x'
	replace sigma_e_n4 = sigma_e`x' if vol_n4==`x'
}

gen signal_2step_n4 = ((sigma_u_n4^2) / ((sigma_u_n4^2) + ((sigma_e_n4^2)/vol)))

gen adj_docfe = docfe*signal_2step_n4 + (1-signal_2step_n4)*mean_docfe

su docfe
return scalar d_std = r(sd)
su adj_docfe
return scalar d_adjstd = r(sd)

su docfe
return scalar d_var = r(Var)
su adj_docfe
return scalar d_adjvar = r(Var)

su docfe, det
return scalar d_p10 = r(p10)
su adj_docfe, det
return scalar d_adjp10 = r(p10)

su docfe, det
return scalar d_p25 = r(p25)
su adj_docfe, det
return scalar d_adjp25 = r(p25)

su docfe, det
return scalar d_p50 = r(p50)
su adj_docfe, det
return scalar d_adjp50 = r(p50)

su docfe, det
return scalar d_p75 = r(p75)
su adj_docfe, det
return scalar d_adjp75 = r(p75)

su docfe, det
return scalar d_p90 = r(p90)
su adj_docfe, det
return scalar d_adjp90 = r(p90)

restore

end


** Then run bootstrap here

#delimit ;
* Simulate 199 times to test;
simulate diff_std = r(d_std) diff_adjstd = r(d_adjstd) diff_var = r(d_var) diff_adjvar = r(d_adjvar)
diff_p10 = r(d_p10) diff_adjp10 = r(d_adjp10) diff_p25 = r(d_p25) diff_adjp25 = r(d_adjp25) diff_p50 = r(d_p50) diff_adjp50 = r(d_adjp50)
diff_p75 = r(d_p75) diff_adjp75 = r(d_adjp75) diff_p90 = r(d_p90) diff_adjp90 = r(d_adjp90)
, reps(199) seed(32786105): myboot2;


#delimit ;
* Boostrap;
bstat, stat(observe_std, observe_adjstd, observe_var, observe_adjvar,
observe_p10, observe_adjp10, observe_p25, observe_adjp25, observe_p50, observe_adjp50,
observe_p75, observe_adjp75, observe_p90, observe_adjp90);

#delimit cr
* Output results
putexcel set "$results/TableA5_col1.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)


* Repeat analysis for 365-day survival (column 2)
**************************************************

use "$savedata/masterdata.dta", replace
keep if sample25==1
gen vol = vol25

cap drop _merge
merge m:1 doctor_id using "$savedata/doctor_vol.dta"

preserve

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

foreach x in 1 2 3 4{
xtreg residual c.std_ami3 c.std_nonami3 if vol_n4==`x', fe i(doctor_id)
predict docfe`x', u
gen sigma_u`x' = e(sigma_u)
gen sigma_e`x' = e(sigma_e)
}

collapse (mean) docfe* residual sigma* vol, by(doctor_id)

xtile vol_n4 = vol, n(4)
egen mean_docfe = mean(docfe), by(vol_n4)

gen sigma_u_n4 = sigma_u1 if vol_n4==1
gen sigma_e_n4 = sigma_e1 if vol_n4==1
foreach x in 2 3 4{
	replace sigma_u_n4 = sigma_u`x' if vol_n4==`x'
	replace sigma_e_n4 = sigma_e`x' if vol_n4==`x'
}

gen signal_2step_n4 = ((sigma_u_n4^2) / ((sigma_u_n4^2) + ((sigma_e_n4^2)/vol)))

gen adj_docfe = docfe*signal_2step_n4 + (1-signal_2step_n4)*mean_docfe

su docfe
matrix observe_std = r(sd)
su adj_docfe
matrix observe_adjstd = r(sd)

su docfe
matrix observe_var = r(Var)
su adj_docfe
matrix observe_adjvar = r(Var)

su docfe, det
matrix observe_p10 = r(p10)
su adj_docfe, det
matrix observe_adjp10 = r(p10)

su docfe, det
matrix observe_p25 = r(p25)
su adj_docfe, det
matrix observe_adjp25 = r(p25)

su docfe, det
matrix observe_p50 = r(p50)
su adj_docfe, det
matrix observe_adjp50 = r(p50)

su docfe, det
matrix observe_p75 = r(p75)
su adj_docfe, det
matrix observe_adjp75 = r(p75)

su docfe, det
matrix observe_p90 = r(p90)
su adj_docfe, det
matrix observe_adjp90 = r(p90)

restore

capture program drop myboot2

program define myboot2, rclass

* draw my bootstrap sample
preserve
bsample, strata(doctor_id)

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

foreach x in 1 2 3 4{
xtreg residual c.std_ami3 c.std_nonami3 if vol_n4==`x', fe i(doctor_id)
predict docfe`x', u
gen sigma_u`x' = e(sigma_u)
gen sigma_e`x' = e(sigma_e)
}

collapse (mean) docfe* residual sigma* vol, by(doctor_id)

xtile vol_n4 = vol, n(4)
egen mean_docfe = mean(docfe), by(vol_n4)

gen sigma_u_n4 = sigma_u1 if vol_n4==1
gen sigma_e_n4 = sigma_e1 if vol_n4==1
foreach x in 2 3 4{
	replace sigma_u_n4 = sigma_u`x' if vol_n4==`x'
	replace sigma_e_n4 = sigma_e`x' if vol_n4==`x'
}

gen signal_2step_n4 = ((sigma_u_n4^2) / ((sigma_u_n4^2) + ((sigma_e_n4^2)/vol)))

gen adj_docfe = docfe*signal_2step_n4 + (1-signal_2step_n4)*mean_docfe

su docfe
return scalar d_std = r(sd)
su adj_docfe
return scalar d_adjstd = r(sd)

su docfe
return scalar d_var = r(Var)
su adj_docfe
return scalar d_adjvar = r(Var)

su docfe, det
return scalar d_p10 = r(p10)
su adj_docfe, det
return scalar d_adjp10 = r(p10)

su docfe, det
return scalar d_p25 = r(p25)
su adj_docfe, det
return scalar d_adjp25 = r(p25)

su docfe, det
return scalar d_p50 = r(p50)
su adj_docfe, det
return scalar d_adjp50 = r(p50)

su docfe, det
return scalar d_p75 = r(p75)
su adj_docfe, det
return scalar d_adjp75 = r(p75)

su docfe, det
return scalar d_p90 = r(p90)
su adj_docfe, det
return scalar d_adjp90 = r(p90)

restore

end


** Then run bootstrap here

#delimit ;
* Simulate 199 times to test;
simulate diff_std = r(d_std) diff_adjstd = r(d_adjstd) diff_var = r(d_var) diff_adjvar = r(d_adjvar)
diff_p10 = r(d_p10) diff_adjp10 = r(d_adjp10) diff_p25 = r(d_p25) diff_adjp25 = r(d_adjp25) diff_p50 = r(d_p50) diff_adjp50 = r(d_adjp50)
diff_p75 = r(d_p75) diff_adjp75 = r(d_adjp75) diff_p90 = r(d_p90) diff_adjp90 = r(d_adjp90)
, reps(199) seed(32786105): myboot2;


#delimit ;
* Boostrap;
bstat, stat(observe_std, observe_adjstd, observe_var, observe_adjvar,
observe_p10, observe_adjp10, observe_p25, observe_adjp25, observe_p50, observe_adjp50,
observe_p75, observe_adjp75, observe_p90, observe_adjp90);

#delimit cr
* Output results
putexcel set "$results/TableA5_col2.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)

